clear all
set more off
set mem 10000000
set matsize 10000
version 15

******************************************************
*** RDROBUST first stage, VD electricity variables ***
******************************************************

** Set file paths
do "$path_code/paths.do"

** Set graph scheme
cd "$path/code/analyze"
set scheme fb, perm

****************************************************************** 
****************************************************************** 

** 1. Main VD_elec RD: RDROBUST, tri kernel, default SE, STFEs, pre lights ctrls 
**                     (1998-2004), drop lights_diff>=20, drop pop_mismatch20, fs
{
use "$panel/panel_dataset_full.dta", clear

	// Keep villages in RD sample
gen in_fs_sample = vplan4<11 & corr_state==1 & sample==1 & sample_h==1 & sing_h==1 & pop_non_zero==1 
keep if in_fs_sample==1
	
	// Create state FEs (since RD robust doesn't let you pass them through)
drop if st_code==32 // Kerala, only 3 villages
tab st_code, gen(STFE)	
drop STFE1 // to avoid collinearity

	// Create district FEs (since RD robust doesn't let you pass them through)
gen stdtFE = stdt
tab stdtFE state     if inlist(stdtFE,62,63,64,70,91,175,180,181,202,308,309,320,490,493,494,495,504,505,508)	
replace stdtFE = 999 if inlist(stdtFE,62,63,64,70,91,175,180,181,202,308,309,320,490,493,494,495,504,505,508)
	// one catch-all district FE for districts with so few in-sample villages that they break rdrobust
gen stdthFE = stdtFE
replace stdthFE = 999 if st_code==29 // Karnataka is almost completely missing hour of commercial power!
tab stdtFE, gen(DTFE)	
drop DTFE1 // to avoid collinearity
tab stdthFE, gen(DTFEh)	
drop DTFEh1 // to avoid collinearity

	// Create block groups (for clustering)
egen stdtbk = group(stdt bk_code)

	// Create lights-difference variable, to identify crazy outliers
gen lights_diff = abs(lights_max2011_hat - lights_max2001_hat)

	// Rename variable to make them easier to loop over
rename vdp_pwr_d_com_11 vd_pwr_d_com_11	
rename vdp_pwr_h_*_avg_11 vd_pwr_h_*_11

	// Create variables to store regresson results
gen fs_step = .
gen yvar = ""
gen ifs = ""
gen control = ""
gen fe = ""
gen kernel = ""
gen bwmethod = ""
gen vce = ""
gen polynomial_order = .
gen beta_conv = .
gen beta_robust = .
gen se_conv = .
gen se_robust = .
gen pval_conv = .
gen pval_robust = .
gen lci_conv = .
gen uci_conv = .
gen lci_robust = .
gen uci_robust = .
gen bw_lo = .
gen bw_hi = .
gen nobs_orig = .
gen nobs_left = .
gen nobs_right = .
gen ndist = . 
gen ymean = .
gen fdesc = ""
local folder = "RDROBUST plots fs vdelec"
local h_wide = 200
local row = 0

	// 1.1 Preferred specification
{
	local fs_step = 1
	local ifs = "in_fs_sample==1 & pop_mismatch20==0 & lights_diff<20"
	local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
	local fe = "STFE*"
	local kernel = "tri"
	local bwmethod = "mserd"
	local vce = ""
	local poly = 1

	foreach y in d_com d_all h_com d_agr d_dom h_agr h_dom h_all {

		// Define outcome variable
		local yvar = "vd_pwr_`y'_11"
		if "`y'"=="d_agr" {
			local title = "Dummy for Agricultural Power Access, 2011"
			local ytitle = "Share of villages, residuals"
			local ygap = ""
		}
		if "`y'"=="d_com" {
			local title = "Dummy for Commercial Power Access, 2011"
			local ytitle = "Share of villages, residuals"
			local ygap = ""
		}
		if "`y'"=="d_dom" {
			local title = "Dummy for Domestic Power Access, 2011"
			local ytitle = "Share of villages, residuals"
			local ygap = ""
		}
		if "`y'"=="d_all" {
			local title = "Dummy for Power Access, All End Uses, 2011"
			local ytitle = "Share of villages, residuals"
			local ygap = ""
		}
		if "`y'"=="h_agr" {
			local title = "Hours of Agricultural Power Access, 2011"
			local ytitle = "Avg hours per day, residuals"
			local ygap = "yscale(titlegap(*-30))"
		}
		if "`y'"=="h_com" {
			local title = "Hours of Commercial Power Access, 2011"
			local ytitle = "Avg hours per day, residuals"
			local ygap = ""
		}
		if "`y'"=="h_dom" {
			local title = "Hours of Domestic Power Access, 2011"
			local ytitle = "Avg hours per day, residuals"
			local ygap = "yscale(titlegap(*-30))"
		}
		if "`y'"=="h_all" {
			local title = "Hours of Power Access, All End Uses, 2011"
			local ytitle = "Avg hours per day, residuals"
			local ygap = ""
		}
		local fdesc = "Preferred"
		local ftag = "preferred"

		// Run first-stage regresssion
		rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

		// Generate in-sample indicator and store bandwidths
		cap drop temp_in_reg
		qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r)) & `yvar'!=.
		local h_l = e(h_l)
		local h_r = e(h_r)

		// Store results
		local row = `row' + 1
		qui replace fs_step = `fs_step' in `row'
		qui replace yvar = "`yvar'" in `row'
		qui replace ifs = "`ifs'" in `row'
		qui replace control = "`controls'" in `row'
		qui replace fe = "`fe'" in `row'
		qui replace kernel = e(kernel) in `row'
		qui replace bwmethod = e(bwselect) in `row'
		qui replace vce = "`vce'" in `row'
		qui replace polynomial_order = e(p) in `row'
		qui replace beta_conv = e(tau_cl) in `row'
		qui replace beta_robust = e(tau_bc) in `row'
		qui replace se_conv = e(se_tau_cl) in `row'
		qui replace se_robust = e(se_tau_rb) in `row'
		qui replace pval_conv = e(pv_cl) in `row'
		qui replace pval_robust = e(pv_rb) in `row'
		qui replace lci_conv = e(ci_l_cl) in `row'
		qui replace uci_conv = e(ci_r_cl) in `row'
		qui replace lci_robust = e(ci_l_rb) in `row'
		qui replace uci_robust = e(ci_r_rb) in `row'
		qui replace bw_lo = e(h_l) in `row'
		qui replace bw_hi = e(h_r) in `row'
		qui replace nobs_orig = e(N) in `row'
		qui replace nobs_left = e(N_h_l) in `row'
		qui replace nobs_right = e(N_h_r) in `row'
		qui unique stdt if temp_in_reg==1
		qui replace ndist = r(unique) in `row'
		qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
		qui replace ymean = r(mean) in `row'
		qui replace fdesc = "`fdesc'" in `row'

		// Residualize `yvar' for regression sample
		qui reg `yvar' `controls' `fe' if temp_in_reg==1
		cap drop y_resid_sample
		predict y_resid_sample, residuals
		
		// Residualize `yvar' for [100,500] bandwidth
		qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
		cap drop y_resid_full
		predict y_resid_full, residuals
				
		// RD plots for regression sample (variable bandwidth) 
		foreach bins in 10 20 40 {
			if `h_l'<=75 {
				local xlab = "250 275 300 325 350"
			}
			else if inrange(`h_l',76,125) {
				local xlab = "200 250 300 350 400"
			}
			else if inrange(`h_l',126,175) {
				local xlab = "150 200 250 300 350 400 450"
			}
			else {
				local xlab = ""
			}
			rdplot y_resid_sample tot_p if temp_in_reg==1, ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`ytitle'", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) `ygap' ///
				xlabel(`xlab', labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/`yvar'_fs_`ftag'_inreg_`bins'.pdf", replace
		}

		// RD plots for regression sample (constant bandwidth) 
		foreach bins in 20 40 {
			rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`ytitle'", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) `ygap' ///
				xlabel(, labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/`yvar'_fs_`ftag'_wide_`bins'.pdf", replace
		}			
	}
}


	// 1.2 Sensitivity to including population mismatches
{		
	local fs_step = 2
	local ifs = "in_fs_sample==1 & lights_diff<20"
	local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
	local fe = "STFE*"
	local kernel = "tri"
	local bwmethod = "mserd"
	local vce = ""
	local poly = 1

	foreach y in d_com d_all h_com h_all {

		// Define outcome variable
		local yvar = "vd_pwr_`y'_11"
		if "`y'"=="d_agr" {
			local title = "2011 Ag Power Access, Pop Mismatches"
			local ytitle = "Share of villages, residuals"
		}
		if "`y'"=="d_com" {
			local title = "2011 Commercial Power Access, Pop Mismatches"
			local ytitle = "Share of villages, residuals"
		}
		if "`y'"=="d_dom" {
			local title = "2011 Domestic Power Access, Pop Mismatches"
			local ytitle = "Share of villages, residuals"
		}
		if "`y'"=="d_all" {
			local title = "2011 All-Use Power Access, Pop Mismatches"
			local ytitle = "Share of villages, residuals"
		}
		if "`y'"=="h_agr" {
			local title = "2011 Ag Power Hours, Pop Mismatches"
			local ytitle = "Avg hours per day, residuals"
		}
		if "`y'"=="h_com" {
			local title = "2011 Commercial Power Hours, Pop Mismatches"
			local ytitle = "Avg hours per day, residuals"
		}
		if "`y'"=="h_dom" {
			local title = "2011 Domestic Power Hours, Pop Mismatches"
			local ytitle = "Avg hours per day, residuals"
		}
		if "`y'"=="h_all" {
			local title = "2011 All-Use Power Hours, Pop Mismatches"
			local ytitle = "Avg hours per day, residuals"
		}
		local fdesc = "popmismatch"
		local ftag = "inclpop20"
		
		// Run first-stage regresssion
		rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

		// Generate in-sample indicator and store bandwidths
		cap drop temp_in_reg
		qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r)) & `yvar'!=.
		local h_l = e(h_l)
		local h_r = e(h_r)
		
		// Store results
		local row = `row' + 1
		qui replace fs_step = `fs_step' in `row'
		qui replace yvar = "`yvar'" in `row'
		qui replace ifs = "`ifs'" in `row'
		qui replace control = "`controls'" in `row'
		qui replace fe = "`fe'" in `row'
		qui replace kernel = e(kernel) in `row'
		qui replace bwmethod = e(bwselect) in `row'
		qui replace vce = "`vce'" in `row'
		qui replace polynomial_order = e(p) in `row'
		qui replace beta_conv = e(tau_cl) in `row'
		qui replace beta_robust = e(tau_bc) in `row'
		qui replace se_conv = e(se_tau_cl) in `row'
		qui replace se_robust = e(se_tau_rb) in `row'
		qui replace pval_conv = e(pv_cl) in `row'
		qui replace pval_robust = e(pv_rb) in `row'
		qui replace lci_conv = e(ci_l_cl) in `row'
		qui replace uci_conv = e(ci_r_cl) in `row'
		qui replace lci_robust = e(ci_l_rb) in `row'
		qui replace uci_robust = e(ci_r_rb) in `row'
		qui replace bw_lo = e(h_l) in `row'
		qui replace bw_hi = e(h_r) in `row'
		qui replace nobs_orig = e(N) in `row'
		qui replace nobs_left = e(N_h_l) in `row'
		qui replace nobs_right = e(N_h_r) in `row'
		qui unique stdt if temp_in_reg==1
		qui replace ndist = r(unique) in `row'
		qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
		qui replace ymean = r(mean) in `row'
		qui replace fdesc = "`fdesc'" in `row'

		// Residualize `yvar' for regression sample
		qui reg `yvar' `controls' `fe' if temp_in_reg==1
		cap drop y_resid_sample
		predict y_resid_sample, residuals
		
		// Residualize `yvar' for [100,500] bandwidth
		qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
		cap drop y_resid_full
		predict y_resid_full, residuals
		
		// RD plots for regression sample (variable bandwidth) 
		foreach bins in 10 20 40 {
			if `h_l'<=75 {
				local xlab = "250 275 300 325 350"
			}
			else if inrange(`h_l',76,125) {
				local xlab = "200 250 300 350 400"
			}
			else if inrange(`h_l',126,175) {
				local xlab = "150 200 250 300 350 400 450"
			}
			else {
				local xlab = ""
			}
			rdplot y_resid_sample tot_p if temp_in_reg==1, ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`ytitle'", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(`xlab', labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/`yvar'_fs_`ftag'_inreg_`bins'.pdf", replace
		}

		// RD plots for regression sample (constant bandwidth) 
		foreach bins in 20 40 {
			rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`ytitle'", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(, labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/`yvar'_fs_`ftag'_wide_`bins'.pdf", replace
		}		
	}
}


	// 1.3 Sensitivity to lights_diff outliers
{	
	local fs_step = 3
	local ifs_base = "in_fs_sample==1 & pop_mismatch20==0"
	local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
	local fe = "STFE*"
	local kernel = "tri"
	local bwmethod = "mserd"
	local vce = ""
	local poly = 1

	foreach y in d_com d_all h_com h_all {

		foreach ld in 10 15 25 30 35 40 45 50  {

			// Define threshould for dropping lights_diff outliers
			local ifs = "`ifs_base' & lights_diff<`ld'"
			local fdesc = "lights outliers `ld'"
			local ftag = "ldiff`ld'"
				
			// Define outcome variable
			local yvar = "vd_pwr_`y'_11"
			if "`y'"=="d_agr" {
				local title = "2011 Ag Power Access, w/o >`ld' Lights Changes"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_com" {
				local title = "2011 Commercial Power Access, w/o >`ld' Lights Changes"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_dom" {
				local title = "2011 Domestic Power Access, w/o >`ld' Lights Changes"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_all" {
				local title = "2011 All-Use Power Access, w/o >`ld' Lights Changes"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="h_agr" {
				local title = "2011 Ag Power Hours, w/o >`ld' Lights Changes"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_com" {
				local title = "2011 Commercial Power Hours, w/o >`ld' Lights Changes"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_dom" {
				local title = "2011 Domestic Power Hours, w/o >`ld' Lights Changes"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_all" {
				local title = "2011 All-Use Power Hours, w/o >`ld' Lights Changes"
				local ytitle = "Avg hours per day, residuals"
			}
			
			// Run first-stage regresssion
			rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

			// Generate in-sample indicator and store bandwidths
			cap drop temp_in_reg
			qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r)) & `yvar'!=.
			local h_l = e(h_l)
			local h_r = e(h_r)
			
			// Store results
			local row = `row' + 1
			qui replace fs_step = `fs_step' in `row'
			qui replace yvar = "`yvar'" in `row'
			qui replace ifs = "`ifs'" in `row'
			qui replace control = "`controls'" in `row'
			qui replace fe = "`fe'" in `row'
			qui replace kernel = e(kernel) in `row'
			qui replace bwmethod = e(bwselect) in `row'
			qui replace vce = "`vce'" in `row'
			qui replace polynomial_order = e(p) in `row'
			qui replace beta_conv = e(tau_cl) in `row'
			qui replace beta_robust = e(tau_bc) in `row'
			qui replace se_conv = e(se_tau_cl) in `row'
			qui replace se_robust = e(se_tau_rb) in `row'
			qui replace pval_conv = e(pv_cl) in `row'
			qui replace pval_robust = e(pv_rb) in `row'
			qui replace lci_conv = e(ci_l_cl) in `row'
			qui replace uci_conv = e(ci_r_cl) in `row'
			qui replace lci_robust = e(ci_l_rb) in `row'
			qui replace uci_robust = e(ci_r_rb) in `row'
			qui replace bw_lo = e(h_l) in `row'
			qui replace bw_hi = e(h_r) in `row'
			qui replace nobs_orig = e(N) in `row'
			qui replace nobs_left = e(N_h_l) in `row'
			qui replace nobs_right = e(N_h_r) in `row'
			qui unique stdt if temp_in_reg==1
			qui replace ndist = r(unique) in `row'
			qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
			qui replace ymean = r(mean) in `row'
			qui replace fdesc = "`fdesc'" in `row'

			// Residualize `yvar' for regression sample
			qui reg `yvar' `controls' `fe' if temp_in_reg==1
			cap drop y_resid_sample
			predict y_resid_sample, residuals
			
			// Residualize `yvar' for [100,500] bandwidth
			qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
			cap drop y_resid_full
			predict y_resid_full, residuals
			
			// RD plots for regression sample (variable bandwidth) 
			foreach bins in 10 20 40 {
				if `h_l'<=75 {
					local xlab = "250 275 300 325 350"
				}
				else if inrange(`h_l',76,125) {
					local xlab = "200 250 300 350 400"
				}
				else if inrange(`h_l',126,175) {
					local xlab = "150 200 250 300 350 400 450"
				}
				else {
					local xlab = ""
				}
				rdplot y_resid_sample tot_p if temp_in_reg==1, ///
					c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
					graph_options( ///
					title("`title'", color(black) size(large)) ///
					ytitle("`ytitle'", size(medlarge)) ///
					xtitle("2001 village population", size(medlarge)) ///
					ylabel(,nogrid angle(0) labsize(medlarge)) ///
					xlabel(`xlab', labsize(medlarge)) ///
					graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
					legend(off))
				graph export "$results/`folder'/`yvar'_fs_`ftag'_inreg_`bins'.pdf", replace
			}

			// RD plots for regression sample (constant bandwidth) 
			foreach bins in 20 40 {
				rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
					c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
					graph_options( ///
					title("`title'", color(black) size(large)) ///
					ytitle("`ytitle'", size(medlarge)) ///
					xtitle("2001 village population", size(medlarge)) ///
					ylabel(,nogrid angle(0) labsize(medlarge)) ///
					xlabel(, labsize(medlarge)) ///
					graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
					legend(off))
				graph export "$results/`folder'/`yvar'_fs_`ftag'_wide_`bins'.pdf", replace
			}		
		}
	}
}


	// 1.4 Sensitivity to different bandwidth calculations
{	
	local fs_step = 4
	local ifs = "in_fs_sample==1 & pop_mismatch20==0 & lights_diff<20"
	local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
	local fe = "STFE*"
	local kernel = "tri"
	local vce = ""
	local poly = 1

	foreach y in d_com d_all h_com h_all {

		foreach reg in 1 2 3 4 5 {

			// Define bandwidth method
			if `reg'==1 {
				local bwmethod = "msetwo"
			}
			if `reg'==2 {
				local bwmethod = "msesum"
			}
			if `reg'==3 {
				local bwmethod = "cerrd"
			}
			if `reg'==4 {
				local bwmethod = "certwo"
			}
			if `reg'==5 {
				local bwmethod = "cersum"
			}
			local fdesc = "BW `bwmethod'"
			local ftag = "BW`bwmethod'"

			// Define outcome variable
			local yvar = "vd_pwr_`y'_11"
			if "`y'"=="d_agr" {
				local title = "2011 Ag Power Access, BW Method `bwmethod'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_com" {
				local title = "2011 Commercial Power Access, BW Method `bwmethod'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_dom" {
				local title = "2011 Domestic Power Access, BW Method `bwmethod'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_all" {
				local title = "2011 All-Use Power Access, BW Method `bwmethod'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="h_agr" {
				local title = "2011 Ag Power Hours, BW Method `bwmethod'"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_com" {
				local title = "2011 Commercial Power Hours, BW Method `bwmethod'"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_dom" {
				local title = "2011 Domestic Power Hours, BW Method `bwmethod'"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_all" {
				local title = "2011 All-Use Power Hours, BW Method `bwmethod'"
				local ytitle = "Avg hours per day, residuals"
			}

			// Run first-stage regresssion
			rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

			// Generate in-sample indicator and store bandwidths
			cap drop temp_in_reg
			qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r)) & `yvar'!=.
			local h_l = e(h_l)
			local h_r = e(h_r)
			
			// Store results
			local row = `row' + 1
			qui replace fs_step = `fs_step' in `row'
			qui replace yvar = "`yvar'" in `row'
			qui replace ifs = "`ifs'" in `row'
			qui replace control = "`controls'" in `row'
			qui replace fe = "`fe'" in `row'
			qui replace kernel = e(kernel) in `row'
			qui replace bwmethod = e(bwselect) in `row'
			qui replace vce = "`vce'" in `row'
			qui replace polynomial_order = e(p) in `row'
			qui replace beta_conv = e(tau_cl) in `row'
			qui replace beta_robust = e(tau_bc) in `row'
			qui replace se_conv = e(se_tau_cl) in `row'
			qui replace se_robust = e(se_tau_rb) in `row'
			qui replace pval_conv = e(pv_cl) in `row'
			qui replace pval_robust = e(pv_rb) in `row'
			qui replace lci_conv = e(ci_l_cl) in `row'
			qui replace uci_conv = e(ci_r_cl) in `row'
			qui replace lci_robust = e(ci_l_rb) in `row'
			qui replace uci_robust = e(ci_r_rb) in `row'
			qui replace bw_lo = e(h_l) in `row'
			qui replace bw_hi = e(h_r) in `row'
			qui replace nobs_orig = e(N) in `row'
			qui replace nobs_left = e(N_h_l) in `row'
			qui replace nobs_right = e(N_h_r) in `row'
			qui unique stdt if temp_in_reg==1
			qui replace ndist = r(unique) in `row'
			qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
			qui replace ymean = r(mean) in `row'
			qui replace fdesc = "`fdesc'" in `row'

			// Residualize `yvar' for regression sample
			qui reg `yvar' `controls' `fe' if temp_in_reg==1
			cap drop y_resid_sample
			predict y_resid_sample, residuals
			
			// RD plots for regression sample (variable bandwidth) 
			foreach bins in 10 20 40 {
				if `h_l'<=75 {
					local xlab = "250 275 300 325 350"
				}
				else if inrange(`h_l',76,125) {
					local xlab = "200 250 300 350 400"
				}
				else if inrange(`h_l',126,175) {
					local xlab = "150 200 250 300 350 400 450"
				}
				else {
					local xlab = ""
				}
				rdplot y_resid_sample tot_p if temp_in_reg==1, ///
					c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
					graph_options( ///
					title("`title'", color(black) size(large)) ///
					ytitle("`ytitle'", size(medlarge)) ///
					xtitle("2001 village population", size(medlarge)) ///
					ylabel(,nogrid angle(0) labsize(medlarge)) ///
					xlabel(`xlab', labsize(medlarge)) ///
					graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
					legend(off))
				graph export "$results/`folder'/`yvar'_fs_`ftag'_inreg_`bins'.pdf", replace
			}
		}		
	}
}


	// 1.5 Sensitivity to 2nd-order polynomial and differet kernels
{	
	local fs_step = 5
	local ifs = "in_fs_sample==1 & pop_mismatch20==0 & lights_diff<20"
	local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
	local fe = "STFE*"
	local bwmethod = "mserd"
	local vce = ""

	foreach y in d_com d_all h_com h_all {

		foreach reg in 1 2 3 4 5 {

			// Define kernel and polynomial combo
			if `reg'==1 {
				local kernel = "tri"
				local poly = 2
				local title_tag = "Tri Kernel, Quadratic"
			}
			if `reg'==2 {
				local kernel = "epa"
				local poly = 1
				local title_tag = "Epa Kernel"
			}
			if `reg'==3 {
				local kernel = "epa"
				local poly = 2
				local title_tag = "Epa Kernel, Quadratic"
			}
			if `reg'==4 {
				local kernel = "uni"
				local poly = 1
				local title_tag = "Uni Kernel"
			}
			if `reg'==5 {
				local kernel = "uni"
				local poly = 2
				local title_tag = "Uni Kernel, Quadratic"
			}
			local fdesc = "kernel `kernel', poly `poly'"
			local ftag = "K`kernel'P`poly'"

			// Define outcome variable
			local yvar = "vd_pwr_`y'_11"
			if "`y'"=="d_agr" {
				local title = "2011 Ag Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_com" {
				local title = "2011 Commercial Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_dom" {
				local title = "2011 Domestic Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_all" {
				local title = "2011 All-Use Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="h_agr" {
				local title = "2011 Ag Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_com" {
				local title = "2011 Commercial Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_dom" {
				local title = "2011 Domestic Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_all" {
				local title = "2011 All-Use Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}

			// Run first-stage regresssion
			rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

			// Generate in-sample indicator and store bandwidths
			cap drop temp_in_reg
			qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r)) & `yvar'!=.
			local h_l = e(h_l)
			local h_r = e(h_r)
			
			// Store results
			local row = `row' + 1
			qui replace fs_step = `fs_step' in `row'
			qui replace yvar = "`yvar'" in `row'
			qui replace ifs = "`ifs'" in `row'
			qui replace control = "`controls'" in `row'
			qui replace fe = "`fe'" in `row'
			qui replace kernel = e(kernel) in `row'
			qui replace bwmethod = e(bwselect) in `row'
			qui replace vce = "`vce'" in `row'
			qui replace polynomial_order = e(p) in `row'
			qui replace beta_conv = e(tau_cl) in `row'
			qui replace beta_robust = e(tau_bc) in `row'
			qui replace se_conv = e(se_tau_cl) in `row'
			qui replace se_robust = e(se_tau_rb) in `row'
			qui replace pval_conv = e(pv_cl) in `row'
			qui replace pval_robust = e(pv_rb) in `row'
			qui replace lci_conv = e(ci_l_cl) in `row'
			qui replace uci_conv = e(ci_r_cl) in `row'
			qui replace lci_robust = e(ci_l_rb) in `row'
			qui replace uci_robust = e(ci_r_rb) in `row'
			qui replace bw_lo = e(h_l) in `row'
			qui replace bw_hi = e(h_r) in `row'
			qui replace nobs_orig = e(N) in `row'
			qui replace nobs_left = e(N_h_l) in `row'
			qui replace nobs_right = e(N_h_r) in `row'
			qui unique stdt if temp_in_reg==1
			qui replace ndist = r(unique) in `row'
			qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
			qui replace ymean = r(mean) in `row'
			qui replace fdesc = "`fdesc'" in `row'

			// Residualize `yvar' for regression sample
			qui reg `yvar' `controls' `fe' if temp_in_reg==1
			cap drop y_resid_sample
			predict y_resid_sample, residuals
			
			// Residualize `yvar' for [100,500] bandwidth
			qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
			cap drop y_resid_full
			predict y_resid_full, residuals
			
			// RD plots for regression sample (variable bandwidth) 
			foreach bins in 10 20 40 {
				if `h_l'<=75 {
					local xlab = "250 275 300 325 350"
				}
				else if inrange(`h_l',76,125) {
					local xlab = "200 250 300 350 400"
				}
				else if inrange(`h_l',126,175) {
					local xlab = "150 200 250 300 350 400 450"
				}
				else {
					local xlab = ""
				}
				rdplot y_resid_sample tot_p if temp_in_reg==1, ///
					c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
					graph_options( ///
					title("`title'", color(black) size(large)) ///
					ytitle("`ytitle'", size(medlarge)) ///
					xtitle("2001 village population", size(medlarge)) ///
					ylabel(,nogrid angle(0) labsize(medlarge)) ///
					xlabel(`xlab', labsize(medlarge)) ///
					graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
					legend(off))
				graph export "$results/`folder'/`yvar'_fs_`ftag'_inreg_`bins'.pdf", replace
			}

			// RD plots for regression sample (constant bandwidth) 
			foreach bins in 20 40 {
				rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
					c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
					graph_options( ///
					title("`title'", color(black) size(large)) ///
					ytitle("`ytitle'", size(medlarge)) ///
					xtitle("2001 village population", size(medlarge)) ///
					ylabel(,nogrid angle(0) labsize(medlarge)) ///
					xlabel(, labsize(medlarge)) ///
					graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
					legend(off))
				graph export "$results/`folder'/`yvar'_fs_`ftag'_wide_`bins'.pdf", replace
			}
		}		
	}
}


	// 1.6 Sensitivity to fixed effects
{	
	local fs_step = 6
	local ifs = "in_fs_sample==1 & pop_mismatch20==0 & lights_diff<20"
	local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
	local kernel = "tri"
	local bwmethod = "mserd"
	local vce = ""
	local poly = 1

	foreach y in d_com d_all h_com h_all {

		foreach reg in 1 2 {

			// Define fixed effects
			if `reg'==1 {
				local fe = ""
				local title_tag = "No FEs"
				local fdesc = "nofes"
				local ftag = "nofes"
			}
			if `reg'==2 & substr("`y'",1,1)=="d" {
				local fe = "DTFE*"
				local title_tag = "District FEs"
				local fdesc = "dtfes"
				local ftag = "dtfes"
			}
			if `reg'==2 & substr("`y'",1,1)=="h" {
				local fe = "DTFEh*"
				local title_tag = "District FEs"
				local fdesc = "dtfes"
				local ftag = "dtfes"
			}

			// Define outcome variable
			local yvar = "vd_pwr_`y'_11"
			if "`y'"=="d_agr" {
				local title = "2011 Ag Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_com" {
				local title = "2011 Commercial Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_dom" {
				local title = "2011 Domestic Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_all" {
				local title = "2011 All-Use Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="h_agr" {
				local title = "2011 Ag Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_com" {
				local title = "2011 Commercial Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_dom" {
				local title = "2011 Domestic Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_all" {
				local title = "2011 All-Use Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}

			// Run first-stage regresssion
			rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

			// Generate in-sample indicator and store bandwidths
			cap drop temp_in_reg
			qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
			local h_l = e(h_l)
			local h_r = e(h_r)

			// Store results
			local row = `row' + 1
			qui replace fs_step = `fs_step' in `row'
			qui replace yvar = "`yvar'" in `row'
			qui replace ifs = "`ifs'" in `row'
			qui replace control = "`controls'" in `row'
			qui replace fe = "`fe'" in `row'
			qui replace kernel = e(kernel) in `row'
			qui replace bwmethod = e(bwselect) in `row'
			qui replace vce = "`vce'" in `row'
			qui replace polynomial_order = e(p) in `row'
			qui replace beta_conv = e(tau_cl) in `row'
			qui replace beta_robust = e(tau_bc) in `row'
			qui replace se_conv = e(se_tau_cl) in `row'
			qui replace se_robust = e(se_tau_rb) in `row'
			qui replace pval_conv = e(pv_cl) in `row'
			qui replace pval_robust = e(pv_rb) in `row'
			qui replace lci_conv = e(ci_l_cl) in `row'
			qui replace uci_conv = e(ci_r_cl) in `row'
			qui replace lci_robust = e(ci_l_rb) in `row'
			qui replace uci_robust = e(ci_r_rb) in `row'
			qui replace bw_lo = e(h_l) in `row'
			qui replace bw_hi = e(h_r) in `row'
			qui replace nobs_orig = e(N) in `row'
			qui replace nobs_left = e(N_h_l) in `row'
			qui replace nobs_right = e(N_h_r) in `row'
			qui unique stdt if temp_in_reg==1
			qui replace ndist = r(unique) in `row'
			qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
			qui replace ymean = r(mean) in `row'
			qui replace fdesc = "`fdesc'" in `row'

			// Residualize `yvar' for regression sample
			qui reg `yvar' `controls' `fe' if temp_in_reg==1
			cap drop y_resid_sample
			predict y_resid_sample, residuals
			
			// Residualize `yvar' for [100,500] bandwidth
			qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
			cap drop y_resid_full
			predict y_resid_full, residuals
			
			// RD plots for regression sample (variable bandwidth) 
			foreach bins in 10 20 40 {
				if `h_l'<=75 {
					local xlab = "250 275 300 325 350"
				}
				else if inrange(`h_l',76,125) {
					local xlab = "200 250 300 350 400"
				}
				else if inrange(`h_l',126,175) {
					local xlab = "150 200 250 300 350 400 450"
				}
				else {
					local xlab = ""
				}
				rdplot y_resid_sample tot_p if temp_in_reg==1, ///
					c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
					graph_options( ///
					title("`title'", color(black) size(large)) ///
					ytitle("`ytitle'", size(medlarge)) ///
					xtitle("2001 village population", size(medlarge)) ///
					ylabel(,nogrid angle(0) labsize(medlarge)) ///
					xlabel(`xlab', labsize(medlarge)) ///
					graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
					legend(off))
				graph export "$results/`folder'/`yvar'_fs_`ftag'_inreg_`bins'.pdf", replace
			}

			// RD plots for regression sample (constant bandwidth) 
			foreach bins in 20 40 {
				rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
					c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
					graph_options( ///
					title("`title'", color(black) size(large)) ///
					ytitle("`ytitle'", size(medlarge)) ///
					xtitle("2001 village population", size(medlarge)) ///
					ylabel(,nogrid angle(0) labsize(medlarge)) ///
					xlabel(, labsize(medlarge)) ///
					graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
					legend(off))
				graph export "$results/`folder'/`yvar'_fs_`ftag'_wide_`bins'.pdf", replace
			}
		}		
	}
}


	// 1.7 Sensitivity to pre controls
{	
	local fs_step = 7
	local ifs = "in_fs_sample==1 & pop_mismatch20==0 & lights_diff<20"
	local fe = "STFE*" 
	local kernel = "tri"
	local bwmethod = "mserd"
	local vce = ""
	local poly = 1

	foreach y in d_com d_all h_com h_all {

		foreach reg in 1 2 3 4 5 6 7 8 /*9 10*/ {

			// Define pre-controls
			if `reg'==1 {
				local controls = "lights_max2001" 
				local title_tag = "2001 Lights Control"
				local fdesc = "ctrl01"
			}
			if `reg'==2 {
				local controls = "lights_max2000 lights_max2001 lights_max2002" 
				local title_tag = "2000-2002 Lights Controls"
				local fdesc = "ctrl0002"
			}
			if `reg'==3 {
				local controls = "lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003" 
				local title_tag = "1999-2003 Lights Controls"
				local fdesc = "ctrl9903"
			}
			if `reg'==4 {
				local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004" 
				local title_tag = "1998-2004 Lights Controls"
				local fdesc = "ctrl9804"
			}
			if `reg'==5 {
				local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005 vdpr_tra_d_app_pr_01" 
				local title_tag = "Control: 2001 road"
				local fdesc = "ctrlroad"
			}
			if `reg'==6 {
				local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005 pct_irr_01" 
				local title_tag = "Control: 2001 share irrigated"
				local fdesc = "ctrlpctirr"
			}
			if `reg'==7 {
				local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005 lit_p_01" 
				local title_tag = "2001 literacy"
				local fdesc = "ctrllit"
			}
			if `reg'==8 {
				local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005 pct_scst_01" 
				local title_tag = "Control: 2001 share SC/ST"
				local fdesc = "ctrlscst"
			}
			if `reg'==9 {
				local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005 vd_pwr_d_any_01" 
				local title_tag = "Control: 2001 power in any sector"
				local fdesc = "ctrlpwrany"
			}
			if `reg'==10 {
				local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005 vd_pwr_d_all_01" 
				local title_tag = "Control: 2001 all-use power"
				local fdesc = "ctrlpwrall"
			}
			local ftag = "`fdesc'"

			// Define outcome variable
			local yvar = "vd_pwr_`y'_11"
			if "`y'"=="d_agr" {
				local title = "2011 Ag Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_com" {
				local title = "2011 Commercial Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_dom" {
				local title = "2011 Domestic Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_all" {
				local title = "2011 All-Use Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="h_agr" {
				local title = "2011 Ag Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_com" {
				local title = "2011 Commercial Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_dom" {
				local title = "2011 Domestic Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_all" {
				local title = "2011 All-Use Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}

			// Run first-stage regresssion
			rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

			// Generate in-sample indicator and store bandwidths
			cap drop temp_in_reg
			qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r)) & `yvar'!=.
			local h_l = e(h_l)
			local h_r = e(h_r)
			
			// Store results
			local row = `row' + 1
			qui replace fs_step = `fs_step' in `row'
			qui replace yvar = "`yvar'" in `row'
			qui replace ifs = "`ifs'" in `row'
			qui replace control = "`controls'" in `row'
			qui replace fe = "`fe'" in `row'
			qui replace kernel = e(kernel) in `row'
			qui replace bwmethod = e(bwselect) in `row'
			qui replace vce = "`vce'" in `row'
			qui replace polynomial_order = e(p) in `row'
			qui replace beta_conv = e(tau_cl) in `row'
			qui replace beta_robust = e(tau_bc) in `row'
			qui replace se_conv = e(se_tau_cl) in `row'
			qui replace se_robust = e(se_tau_rb) in `row'
			qui replace pval_conv = e(pv_cl) in `row'
			qui replace pval_robust = e(pv_rb) in `row'
			qui replace lci_conv = e(ci_l_cl) in `row'
			qui replace uci_conv = e(ci_r_cl) in `row'
			qui replace lci_robust = e(ci_l_rb) in `row'
			qui replace uci_robust = e(ci_r_rb) in `row'
			qui replace bw_lo = e(h_l) in `row'
			qui replace bw_hi = e(h_r) in `row'
			qui replace nobs_orig = e(N) in `row'
			qui replace nobs_left = e(N_h_l) in `row'
			qui replace nobs_right = e(N_h_r) in `row'
			qui unique stdt if temp_in_reg==1
			qui replace ndist = r(unique) in `row'
			qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
			qui replace ymean = r(mean) in `row'
			qui replace fdesc = "`fdesc'" in `row'

			// Residualize `yvar' for regression sample
			qui reg `yvar' `controls' `fe' if temp_in_reg==1
			cap drop y_resid_sample
			predict y_resid_sample, residuals
			
			// Residualize `yvar' for [100,500] bandwidth
			qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
			cap drop y_resid_full
			predict y_resid_full, residuals
			
			// RD plots for regression sample (variable bandwidth) 
			foreach bins in 10 20 40 {
				if `h_l'<=75 {
					local xlab = "250 275 300 325 350"
				}
				else if inrange(`h_l',76,125) {
					local xlab = "200 250 300 350 400"
				}
				else if inrange(`h_l',126,175) {
					local xlab = "150 200 250 300 350 400 450"
				}
				else {
					local xlab = ""
				}
				rdplot y_resid_sample tot_p if temp_in_reg==1, ///
					c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
					graph_options( ///
					title("`title'", color(black) size(large)) ///
					ytitle("`ytitle'", size(medlarge)) ///
					xtitle("2001 village population", size(medlarge)) ///
					ylabel(,nogrid angle(0) labsize(medlarge)) ///
					xlabel(`xlab', labsize(medlarge)) ///
					graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
					legend(off))
				graph export "$results/`folder'/`yvar'_fs_`ftag'_inreg_`bins'.pdf", replace
			}

			// RD plots for regression sample (constant bandwidth) 
			foreach bins in 20 40 {
				rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
					c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
					graph_options( ///
					title("`title'", color(black) size(large)) ///
					ytitle("`ytitle'", size(medlarge)) ///
					xtitle("2001 village population", size(medlarge)) ///
					ylabel(,nogrid angle(0) labsize(medlarge)) ///
					xlabel(, labsize(medlarge)) ///
					graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
					legend(off))
				graph export "$results/`folder'/`yvar'_fs_`ftag'_wide_`bins'.pdf", replace
			}
		}		
	}
}


	// 1.8 Sensitivity to alternative standard errors
{	
	local fs_step = 8
	local ifs = "in_fs_sample==1 & pop_mismatch20==0 & lights_diff<20"
	local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
	local fe = "STFE*"
	local kernel = "tri"
	local bwmethod = "mserd"
	local poly = 1

	foreach y in d_com d_all h_com h_all {

		foreach reg in 1 2 3 4 5 6 7 {

			// Define vce method
			if `reg'==1 {
				local vce = "vce(nn 5)"
				local title_tag = "NN SEs (min 5)"
				local fdesc = "VCEnn5"
			}
			if `reg'==2 {
				local vce = "vce(nncluster stdt)"
				local title_tag = "NN-Clustered SEs (District)"
				local fdesc = "VCEnnclDT"
			}
			if `reg'==3 {
				local vce = "vce(nncluster stdtbk)"
				local title_tag = "NN-Clustered SEs (Block)"
				local fdesc = "VCEnnclBK"
			}
			if `reg'==4 {
				local vce = "vce(nncluster stdt 5)"
				local title_tag = "NN-Clustered SEs (District, min 5)"
				local fdesc = "VCEnnclDT5"
			}
			if `reg'==5 {
				local vce = "vce(nncluster stdtbk 5)"
				local title_tag = "NN-Clustered SEs (Block, min 5)"
				local fdesc = "VCEnnclBK5"
			}
			if `reg'==6 {
				local vce = "vce(cluster stdt)"
				local title_tag = "Clustered SEs (District)"
				local fdesc = "VCEclDT"
			}
			if `reg'==7 {
				local vce = "vce(cluster stdtbk)"
				local title_tag = "Clustered SEs (Block)"
				local fdesc = "VCEclBK"
			}
			local ftag = "`fdesc'"

			// Define outcome variable
			local yvar = "vd_pwr_`y'_11"
			if "`y'"=="d_agr" {
				local title = "2011 Ag Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_com" {
				local title = "2011 Commercial Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_dom" {
				local title = "2011 Domestic Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="d_all" {
				local title = "2011 All-Use Power Access, `title_tag'"
				local ytitle = "Share of villages, residuals"
			}
			if "`y'"=="h_agr" {
				local title = "2011 Ag Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_com" {
				local title = "2011 Commercial Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_dom" {
				local title = "2011 Domestic Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}
			if "`y'"=="h_all" {
				local title = "2011 All-Use Power Hours, `title_tag'"
				local ytitle = "Avg hours per day, residuals"
			}

			// Run first-stage regresssion
			rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

			// Generate in-sample indicator and store bandwidths
			cap drop temp_in_reg
			qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r)) & `yvar'!=.
			local h_l = e(h_l)
			local h_r = e(h_r)
			
			// Store results
			local row = `row' + 1
			qui replace fs_step = `fs_step' in `row'
			qui replace yvar = "`yvar'" in `row'
			qui replace ifs = "`ifs'" in `row'
			qui replace control = "`controls'" in `row'
			qui replace fe = "`fe'" in `row'
			qui replace kernel = e(kernel) in `row'
			qui replace bwmethod = e(bwselect) in `row'
			qui replace vce = "`vce'" in `row'
			qui replace polynomial_order = e(p) in `row'
			qui replace beta_conv = e(tau_cl) in `row'
			qui replace beta_robust = e(tau_bc) in `row'
			qui replace se_conv = e(se_tau_cl) in `row'
			qui replace se_robust = e(se_tau_rb) in `row'
			qui replace pval_conv = e(pv_cl) in `row'
			qui replace pval_robust = e(pv_rb) in `row'
			qui replace lci_conv = e(ci_l_cl) in `row'
			qui replace uci_conv = e(ci_r_cl) in `row'
			qui replace lci_robust = e(ci_l_rb) in `row'
			qui replace uci_robust = e(ci_r_rb) in `row'
			qui replace bw_lo = e(h_l) in `row'
			qui replace bw_hi = e(h_r) in `row'
			qui replace nobs_orig = e(N) in `row'
			qui replace nobs_left = e(N_h_l) in `row'
			qui replace nobs_right = e(N_h_r) in `row'
			qui unique stdt if temp_in_reg==1
			qui replace ndist = r(unique) in `row'
			qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
			qui replace ymean = r(mean) in `row'
			qui replace fdesc = "`fdesc'" in `row'

		}		
	}
}


	// Save results
keep fs_step-fdesc
dropmiss, obs force
compress
save "$results/RDROBUST_fs_vdelec.dta", replace

}

****************************************************************** 
****************************************************************** 

** 2. Falsification tests for main VD_elec RD
{

	// 2.1 11th Plan, single habitation
{
use "$panel/panel_dataset_full.dta", clear

	// Keep villages in FALSIFICATION sample: 11th Plan, single habitation
gen in_fs_sampleFALS = vplan4==11 & corr_state==1 & sample==1 & sample_h==1 & sing_h==1 & pop_non_zero==1 
keep if in_fs_sampleFALS==1
	
	// Create state FEs (since RD robust doesn't let you pass them through)
drop if st_code==32 // Kerala, only 3 villages
tab st_code, gen(STFE)	
drop STFE1 // to avoid collinearity
gen st_codeH = st_code
tab st_codeH, gen(STFEH)	
drop STFEH1 STFEH13

	// Create block groups (for clustering)
egen stdtbk = group(stdt bk_code)

	// Create lights-difference variable, to identify crazy outliers
gen lights_diff = abs(lights_max2011_hat - lights_max2001_hat)

	// Rename variable to make them easier to loop over
rename vdp_pwr_d_com_11 vd_pwr_d_com_11	
rename vdp_pwr_h_*_avg_11 vd_pwr_h_*_11

	// Create variables to store regresson results
gen fs_step = .
gen yvar = ""
gen ifs = ""
gen control = ""
gen fe = ""
gen kernel = ""
gen bwmethod = ""
gen vce = ""
gen polynomial_order = .
gen beta_conv = .
gen beta_robust = .
gen se_conv = .
gen se_robust = .
gen pval_conv = .
gen pval_robust = .
gen lci_conv = .
gen uci_conv = .
gen lci_robust = .
gen uci_robust = .
gen bw_lo = .
gen bw_hi = .
gen nobs_orig = .
gen nobs_left = .
gen nobs_right = .
gen ndist = . 
gen ymean = .
gen fdesc = ""
gen fals = ""
local folder = "RDROBUST plots fs vdelec"
local h_wide = 200

local fs_step = 9
local ifs = "in_fs_sampleFALS==1 & pop_mismatch20==0 & lights_diff<20"
local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "mserd"
local vce = ""
local poly = 1
local fals = "11th Plan, single hab"
local ftag = "FALS_11th_singH"
local row = 0

foreach y in d_com d_all h_com h_all {

	// Define outcome variable
	local yvar = "vd_pwr_`y'_11"
	if "`y'"=="d_com" {
		local title = "2011 Commercial Power Access, 11th Plan, Single Habitation"
		local ytitle = "Share of villages, residuals"
	}
	if "`y'"=="d_all" {
		local title = "2011 All-Uses Power Access, 11th Plan, Single Habitation"
		local ytitle = "Share of villages, residuals"
	}
	if "`y'"=="h_com" {
		local title = "2011 Commercial Hours of Power, 11th Plan, Single Habitation"
		local ytitle = "Avg hours per day, residuals"
		local fe = subinstr("`fe'","FE","FEH",1)
	}

	// Run first-stage regresssion
	rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

	// Generate in-sample indicator and store bandwidths
	cap drop temp_in_reg
	qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r)) & `yvar'!=.
	local h_l = e(h_l)
	local h_r = e(h_r)

	// Store results
	local row = `row' + 1
	qui replace fs_step = `fs_step' in `row'
	qui replace yvar = "`yvar'" in `row'
	qui replace ifs = "`ifs'" in `row'
	qui replace control = "`controls'" in `row'
	qui replace fe = "`fe'" in `row'
	qui replace kernel = e(kernel) in `row'
	qui replace bwmethod = e(bwselect) in `row'
	qui replace vce = "`vce'" in `row'
	qui replace polynomial_order = e(p) in `row'
	qui replace beta_conv = e(tau_cl) in `row'
	qui replace beta_robust = e(tau_bc) in `row'
	qui replace se_conv = e(se_tau_cl) in `row'
	qui replace se_robust = e(se_tau_rb) in `row'
	qui replace pval_conv = e(pv_cl) in `row'
	qui replace pval_robust = e(pv_rb) in `row'
	qui replace lci_conv = e(ci_l_cl) in `row'
	qui replace uci_conv = e(ci_r_cl) in `row'
	qui replace lci_robust = e(ci_l_rb) in `row'
	qui replace uci_robust = e(ci_r_rb) in `row'
	qui replace bw_lo = e(h_l) in `row'
	qui replace bw_hi = e(h_r) in `row'
	qui replace nobs_orig = e(N) in `row'
	qui replace nobs_left = e(N_h_l) in `row'
	qui replace nobs_right = e(N_h_r) in `row'
	qui unique stdt if temp_in_reg==1
	qui replace ndist = r(unique) in `row'
	qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
	qui replace ymean = r(mean) in `row'
	qui replace fdesc = "`ftag'" in `row'
	qui replace fals = "`fals'" in `row'

	// Residualize `yvar' for regression sample
	qui reg `yvar' `controls' `fe' if temp_in_reg==1
	cap drop y_resid_sample
	predict y_resid_sample, residuals
	
	// Residualize `yvar' for [100,500] bandwidth
	qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
	cap drop y_resid_full
	predict y_resid_full, residuals
	
	// RD plots for regression sample (variable bandwidth) 
	foreach bins in 10 20 40 {
		if `h_l'<=75 {
			local xlab = "250 275 300 325 350"
		}
		else if inrange(`h_l',76,125) {
			local xlab = "200 250 300 350 400"
		}
		else if inrange(`h_l',126,175) {
			local xlab = "150 200 250 300 350 400 450"
		}
		else {
			local xlab = ""
		}
		rdplot y_resid_sample tot_p if temp_in_reg==1, ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`ytitle'", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(,nogrid angle(0) labsize(medlarge)) ///
			xlabel(`xlab', labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/`yvar'_fs_`ftag'_inreg_`bins'.pdf", replace
	}

	// RD plots for regression sample (constant bandwidth) 
	foreach bins in 20 40 {
		rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`ytitle'", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(,nogrid angle(0) labsize(medlarge)) ///
			xlabel(, labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/`yvar'_fs_`ftag'_wide_`bins'.pdf", replace
	}			
}

	// Store results in temp file
keep fs_step-fdesc fals
dropmiss, obs force
tempfile fals1
save `fals1'
}


	// 2.2 10th Plan, multi habitation
{
use "$panel/panel_dataset_full.dta", clear

	// Keep villages in FALSIFICATION sample: 10th Plan, multi habitation
gen in_fs_sampleFALS = vplan4<11 & corr_state==1 & sample==1 & sample_h==1 & sing_h==0 & pop_non_zero==1 
keep if in_fs_sampleFALS==1
	
	// Create state FEs (since RD robust doesn't let you pass them through)
drop if st_code==32 // Kerala, only 3 villages
tab st_code, gen(STFE)	
drop STFE1 STFE2 // to avoid collinearity
gen st_codeH = st_code
tab st_codeH, gen(STFEH)	
drop STFEH1 STFEH2 STFEH12

	// Create block groups (for clustering)
egen stdtbk = group(stdt bk_code)

	// Create lights-difference variable, to identify crazy outliers
gen lights_diff = abs(lights_max2011_hat - lights_max2001_hat)

	// Rename variable to make them easier to loop over
rename vdp_pwr_d_com_11 vd_pwr_d_com_11	
rename vdp_pwr_h_*_avg_11 vd_pwr_h_*_11

	// Create variables to store regresson results
gen fs_step = .
gen yvar = ""
gen ifs = ""
gen control = ""
gen fe = ""
gen kernel = ""
gen bwmethod = ""
gen vce = ""
gen polynomial_order = .
gen beta_conv = .
gen beta_robust = .
gen se_conv = .
gen se_robust = .
gen pval_conv = .
gen pval_robust = .
gen lci_conv = .
gen uci_conv = .
gen lci_robust = .
gen uci_robust = .
gen bw_lo = .
gen bw_hi = .
gen nobs_orig = .
gen nobs_left = .
gen nobs_right = .
gen ndist = . 
gen ymean = .
gen fdesc = ""
gen fals = ""
local folder = "RDROBUST plots fs vdelec"
local h_wide = 200

local fs_step = 9
local ifs = "in_fs_sampleFALS==1 & pop_mismatch20==0 & lights_diff<20"
local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "mserd"
local vce = ""
local poly = 1
local fals = "10th Plan, multi hab"
local ftag = "FALS_10th_multiH"
local row = 0

foreach y in d_com d_all h_com h_all {

	// Define outcome variable
	local yvar = "vd_pwr_`y'_11"
	if "`y'"=="d_com" {
		local title = "2011 Commercial Power Access, 10th Plan, Multiple Habitations"
		local ytitle = "Share of villages, residuals"
	}
	if "`y'"=="d_all" {
		local title = "2011 All-Uses Power Access, 10th Plan, Multiple Habitations"
		local ytitle = "Share of villages, residuals"
	}
	if "`y'"=="h_com" {
		local title = "2011 Commercial Hours of Power, 10th Plan, Multiple Habitations"
		local ytitle = "Avg hours per day, residuals"
		local fe = subinstr("`fe'","FE","FEH",1)
	}

	// Run first-stage regresssion
	rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

	// Generate in-sample indicator and store bandwidths
	cap drop temp_in_reg
	qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r)) & `yvar'!=.
	local h_l = e(h_l)
	local h_r = e(h_r)

	// Store results
	local row = `row' + 1
	qui replace fs_step = `fs_step' in `row'
	qui replace yvar = "`yvar'" in `row'
	qui replace ifs = "`ifs'" in `row'
	qui replace control = "`controls'" in `row'
	qui replace fe = "`fe'" in `row'
	qui replace kernel = e(kernel) in `row'
	qui replace bwmethod = e(bwselect) in `row'
	qui replace vce = "`vce'" in `row'
	qui replace polynomial_order = e(p) in `row'
	qui replace beta_conv = e(tau_cl) in `row'
	qui replace beta_robust = e(tau_bc) in `row'
	qui replace se_conv = e(se_tau_cl) in `row'
	qui replace se_robust = e(se_tau_rb) in `row'
	qui replace pval_conv = e(pv_cl) in `row'
	qui replace pval_robust = e(pv_rb) in `row'
	qui replace lci_conv = e(ci_l_cl) in `row'
	qui replace uci_conv = e(ci_r_cl) in `row'
	qui replace lci_robust = e(ci_l_rb) in `row'
	qui replace uci_robust = e(ci_r_rb) in `row'
	qui replace bw_lo = e(h_l) in `row'
	qui replace bw_hi = e(h_r) in `row'
	qui replace nobs_orig = e(N) in `row'
	qui replace nobs_left = e(N_h_l) in `row'
	qui replace nobs_right = e(N_h_r) in `row'
	qui unique stdt if temp_in_reg==1
	qui replace ndist = r(unique) in `row'
	qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
	qui replace ymean = r(mean) in `row'
	qui replace fdesc = "`ftag'" in `row'
	qui replace fals = "`fals'" in `row'

	// Residualize `yvar' for regression sample
	qui reg `yvar' `controls' `fe' if temp_in_reg==1
	cap drop y_resid_sample
	predict y_resid_sample, residuals
	
	// Residualize `yvar' for [100,500] bandwidth
	qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
	cap drop y_resid_full
	predict y_resid_full, residuals
	
	// RD plots for regression sample (variable bandwidth) 
	foreach bins in 10 20 40 {
		if `h_l'<=75 {
			local xlab = "250 275 300 325 350"
		}
		else if inrange(`h_l',76,125) {
			local xlab = "200 250 300 350 400"
		}
		else if inrange(`h_l',126,175) {
			local xlab = "150 200 250 300 350 400 450"
		}
		else {
			local xlab = ""
		}
		rdplot y_resid_sample tot_p if temp_in_reg==1, ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`ytitle'", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(,nogrid angle(0) labsize(medlarge)) ///
			xlabel(`xlab', labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/`yvar'_fs_`ftag'_inreg_`bins'.pdf", replace
	}

	// RD plots for regression sample (constant bandwidth) 
	foreach bins in 20 40 {
		rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`ytitle'", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(,nogrid angle(0) labsize(medlarge)) ///
			xlabel(, labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/`yvar'_fs_`ftag'_wide_`bins'.pdf", replace
	}			
}

	// Store results in temp file
keep fs_step-fdesc fals
dropmiss, obs force
tempfile fals2
save `fals2'
}


	// 2.3 11th Plan, multi habitation
{
use "$panel/panel_dataset_full.dta", clear

	// Keep villages in FALSIFICATION sample: 11th Plan, multi habitation
gen in_fs_sampleFALS = vplan4==11 & corr_state==1 & sample==1 & sample_h==1 & sing_h==0 & pop_non_zero==1 
keep if in_fs_sampleFALS==1
	
	// Create state FEs (since RD robust doesn't let you pass them through)
drop if st_code==32 // Kerala, only 3 villages
tab st_code, gen(STFE)	
drop STFE1 STFE2 STFE4 // to avoid collinearity
gen st_codeH = st_code
tab st_codeH, gen(STFEH)	
drop STFEH1 STFEH2 STFEH4 STFEH13

	// Create block groups (for clustering)
egen stdtbk = group(stdt bk_code)

	// Create lights-difference variable, to identify crazy outliers
gen lights_diff = abs(lights_max2011_hat - lights_max2001_hat)

	// Rename variable to make them easier to loop over
rename vdp_pwr_d_com_11 vd_pwr_d_com_11	
rename vdp_pwr_h_*_avg_11 vd_pwr_h_*_11

	// Create variables to store regresson results
gen fs_step = .
gen yvar = ""
gen ifs = ""
gen control = ""
gen fe = ""
gen kernel = ""
gen bwmethod = ""
gen vce = ""
gen polynomial_order = .
gen beta_conv = .
gen beta_robust = .
gen se_conv = .
gen se_robust = .
gen pval_conv = .
gen pval_robust = .
gen lci_conv = .
gen uci_conv = .
gen lci_robust = .
gen uci_robust = .
gen bw_lo = .
gen bw_hi = .
gen nobs_orig = .
gen nobs_left = .
gen nobs_right = .
gen ndist = . 
gen ymean = .
gen fdesc = ""
gen fals = ""
local folder = "RDROBUST plots fs vdelec"
local h_wide = 200

local fs_step = 9
local ifs = "in_fs_sampleFALS==1 & pop_mismatch20==0 & lights_diff<20"
local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "mserd"
local vce = ""
local poly = 1
local fals = "11th Plan, multi hab"
local ftag = "FALS_11th_multiH"
local row = 0

foreach y in d_com d_all h_com h_all {

	// Define outcome variable
	local yvar = "vd_pwr_`y'_11"
	if "`y'"=="d_com" {
		local title = "2011 Commercial Power Access, 11th Plan, Multiple Habitations"
		local ytitle = "Share of villages, residuals"
	}
	if "`y'"=="d_all" {
		local title = "2011 All-Uses Power Access, 11th Plan, Multiple Habitations"
		local ytitle = "Share of villages, residuals"
	}
	if "`y'"=="h_com" {
		local title = "2011 Commercial Hours of Power, 11th Plan, Multiple Habitations"
		local ytitle = "Avg hours per day, residuals"
		local fe = subinstr("`fe'","FE","FEH",1)
	}

	// Run first-stage regresssion
	rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

	// Generate in-sample indicator and store bandwidths
	cap drop temp_in_reg
	qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r)) & `yvar'!=.
	local h_l = e(h_l)
	local h_r = e(h_r)

	// Store results
	local row = `row' + 1
	qui replace fs_step = `fs_step' in `row'
	qui replace yvar = "`yvar'" in `row'
	qui replace ifs = "`ifs'" in `row'
	qui replace control = "`controls'" in `row'
	qui replace fe = "`fe'" in `row'
	qui replace kernel = e(kernel) in `row'
	qui replace bwmethod = e(bwselect) in `row'
	qui replace vce = "`vce'" in `row'
	qui replace polynomial_order = e(p) in `row'
	qui replace beta_conv = e(tau_cl) in `row'
	qui replace beta_robust = e(tau_bc) in `row'
	qui replace se_conv = e(se_tau_cl) in `row'
	qui replace se_robust = e(se_tau_rb) in `row'
	qui replace pval_conv = e(pv_cl) in `row'
	qui replace pval_robust = e(pv_rb) in `row'
	qui replace lci_conv = e(ci_l_cl) in `row'
	qui replace uci_conv = e(ci_r_cl) in `row'
	qui replace lci_robust = e(ci_l_rb) in `row'
	qui replace uci_robust = e(ci_r_rb) in `row'
	qui replace bw_lo = e(h_l) in `row'
	qui replace bw_hi = e(h_r) in `row'
	qui replace nobs_orig = e(N) in `row'
	qui replace nobs_left = e(N_h_l) in `row'
	qui replace nobs_right = e(N_h_r) in `row'
	qui unique stdt if temp_in_reg==1
	qui replace ndist = r(unique) in `row'
	qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
	qui replace ymean = r(mean) in `row'
	qui replace fdesc = "`ftag'" in `row'
	qui replace fals = "`fals'" in `row'

	// Residualize `yvar' for regression sample
	qui reg `yvar' `controls' `fe' if temp_in_reg==1
	cap drop y_resid_sample
	predict y_resid_sample, residuals
	
	// Residualize `yvar' for [100,500] bandwidth
	qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
	cap drop y_resid_full
	predict y_resid_full, residuals
	
	// RD plots for regression sample (variable bandwidth) 
	foreach bins in 10 20 40 {
		if `h_l'<=75 {
			local xlab = "250 275 300 325 350"
		}
		else if inrange(`h_l',76,125) {
			local xlab = "200 250 300 350 400"
		}
		else if inrange(`h_l',126,175) {
			local xlab = "150 200 250 300 350 400 450"
		}
		else {
			local xlab = ""
		}
		rdplot y_resid_sample tot_p if temp_in_reg==1, ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`ytitle'", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(,nogrid angle(0) labsize(medlarge)) ///
			xlabel(`xlab', labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/`yvar'_fs_`ftag'_inreg_`bins'.pdf", replace
	}

	// RD plots for regression sample (constant bandwidth) 
	foreach bins in 20 40 {
		rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`ytitle'", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(,nogrid angle(0) labsize(medlarge)) ///
			xlabel(, labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/`yvar'_fs_`ftag'_wide_`bins'.pdf", replace
	}			
}

	// Store results in temp file
keep fs_step-fdesc fals
dropmiss, obs force
tempfile fals3
save `fals3'
}


	// Append falsification results to VD elec first stage results file
use "$results/RDROBUST_fs_vdelec.dta", replace
append using `fals1' `fals2' `fals3' 
duplicates drop
compress
save "$results/RDROBUST_fs_vdelec.dta", replace

}

****************************************************************** 
****************************************************************** 

** 3. Adhoc bandwidth sensitivity analysis (manual bandwidths)
{
use "$panel/panel_dataset_full.dta", clear

	// Keep villages in RD sample
gen in_fs_sample = vplan4<11 & corr_state==1 & sample==1 & sample_h==1 & sing_h==1 & pop_non_zero==1 
keep if in_fs_sample==1
	
	// Create state FEs (since RD robust doesn't let you pass them through)
drop if st_code==32 // Kerala, only 3 villages
tab st_code, gen(STFE)	
drop STFE1 // to avoid collinearity

	// Create district FEs (since RD robust doesn't let you pass them through)
gen stdtFE = stdt
tab stdtFE state     if inlist(stdtFE,62,63,64,70,91,175,180,181,202,308,309,320,490,493,494,495,504,505,508)	
replace stdtFE = 999 if inlist(stdtFE,62,63,64,70,91,175,180,181,202,308,309,320,490,493,494,495,504,505,508)
	// one catch-all district FE for districts with so few in-sample villages that they break rdrobust
gen stdthFE = stdtFE
replace stdthFE = 999 if st_code==29 // Karnataka is almost completely missing hour of commercial power!
tab stdtFE, gen(DTFE)	
drop DTFE1 // to avoid collinearity
tab stdthFE, gen(DTFEh)	
drop DTFEh1 // to avoid collinearity

	// Create block groups (for clustering)
egen stdtbk = group(stdt bk_code)

	// Create lights-difference variable, to identify crazy outliers
gen lights_diff = abs(lights_max2011_hat - lights_max2001_hat)

	// Rename variable to make them easier to loop over
rename vdp_pwr_d_com_11 vd_pwr_d_com_11	
rename vdp_pwr_h_*_avg_11 vd_pwr_h_*_11

	// Create variables to store regresson results
gen fs_step = .
gen yvar = ""
gen ifs = ""
gen control = ""
gen fe = ""
gen kernel = ""
gen bwmethod = ""
gen vce = ""
gen polynomial_order = .
gen beta_conv = .
gen beta_robust = .
gen se_conv = .
gen se_robust = .
gen pval_conv = .
gen pval_robust = .
gen lci_conv = .
gen uci_conv = .
gen lci_robust = .
gen uci_robust = .
gen bw_lo = .
gen bw_hi = .
gen nobs_orig = .
gen nobs_left = .
gen nobs_right = .
gen ndist = . 
gen ymean = .
gen fdesc = ""
local folder = "RDROBUST plots fs vdelec"
local row = 0

	// 3.1 JHard-coded bandwidths
{
	local fs_step = 10
	local ifs = "in_fs_sample==1 & pop_mismatch20==0 & lights_diff<20"
	local controls = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
	local fe = "STFE*"
	local kernel = "tri"
	local bwmethod = "manual"
	local vce = ""
	local poly = 1

	foreach y in d_com d_all h_com d_agr d_dom h_agr h_dom h_all {

		// Define outcome variable
		local yvar = "vd_pwr_`y'_11"
		if "`y'"=="d_agr" {
			local title = "Dummy for Agricultural Power Access, 2011"
			local ytitle = "Share of villages, residuals"
			local ygap = "yscale(titlegap(*-30))"
		}
		if "`y'"=="d_com" {
			local title = "Dummy for Commercial Power Access, 2011"
			local ytitle = "Share of villages, residuals"
			local ygap = ""
		}
		if "`y'"=="d_dom" {
			local title = "Dummy for Domestic Power Access, 2011"
			local ytitle = "Share of villages, residuals"
			local ygap = "yscale(titlegap(*-30))"
		}
		if "`y'"=="d_all" {
			local title = "Dummy for Power Access, All End Uses, 2011"
			local ytitle = "Share of villages, residuals"
			local ygap = ""
		}
		if "`y'"=="h_agr" {
			local title = "Hours of Agricultural Power Access, 2011"
			local ytitle = "Avg hours per day, residuals"
			local ygap = ""
		}
		if "`y'"=="h_com" {
			local title = "Hours of Commercial Power Access, 2011"
			local ytitle = "Avg hours per day, residuals"
			local ygap = ""
		}
		if "`y'"=="h_dom" {
			local title = "Hours of Domestic Power Access, 2011"
			local ytitle = "Avg hours per day, residuals"
			local ygap = ""
		}
		if "`y'"=="h_all" {
			local title = "Hours of Power Access, All End Uses, 2011"
			local ytitle = "Avg hours per day, residuals"
			local ygap = ""
		}
		local fdesc = "Manual bandwidth"

		foreach h in 50 60 70 80 90 100 110 120 130 140 150 {
		    
			// Run first-stage regresssion
			rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') h(`h' `h') p(`poly') all `vce'

			// Generate in-sample indicator and store bandwidths
			cap drop temp_in_reg
			qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r)) & `yvar'!=.
			local h_l = e(h_l)
			local h_r = e(h_r)

			// Store results
			local row = `row' + 1
			qui replace fs_step = `fs_step' in `row'
			qui replace yvar = "`yvar'" in `row'
			qui replace ifs = "`ifs'" in `row'
			qui replace control = "`controls'" in `row'
			qui replace fe = "`fe'" in `row'
			qui replace kernel = e(kernel) in `row'
			qui replace bwmethod = "manual" in `row'
			qui replace vce = "`vce'" in `row'
			qui replace polynomial_order = e(p) in `row'
			qui replace beta_conv = e(tau_cl) in `row'
			qui replace beta_robust = e(tau_bc) in `row'
			qui replace se_conv = e(se_tau_cl) in `row'
			qui replace se_robust = e(se_tau_rb) in `row'
			qui replace pval_conv = e(pv_cl) in `row'
			qui replace pval_robust = e(pv_rb) in `row'
			qui replace lci_conv = e(ci_l_cl) in `row'
			qui replace uci_conv = e(ci_r_cl) in `row'
			qui replace lci_robust = e(ci_l_rb) in `row'
			qui replace uci_robust = e(ci_r_rb) in `row'
			qui replace bw_lo = e(h_l) in `row'
			qui replace bw_hi = e(h_r) in `row'
			qui replace nobs_orig = e(N) in `row'
			qui replace nobs_left = e(N_h_l) in `row'
			qui replace nobs_right = e(N_h_r) in `row'
			qui unique stdt if temp_in_reg==1
			qui replace ndist = r(unique) in `row'
			qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
			qui replace ymean = r(mean) in `row'
			qui replace fdesc = "`fdesc'" in `row'

			// Residualize `yvar' for regression sample
			qui reg `yvar' `controls' `fe' if temp_in_reg==1
			cap drop y_resid_sample
			predict y_resid_sample, residuals
								
			// RD plots for regression sample (variable bandwidth) 
			foreach bins in 10 20 40 {
				if `h_l'<=75 {
					local xlab = "250 275 300 325 350"
				}
				else if inrange(`h_l',76,125) {
					local xlab = "200 250 300 350 400"
				}
				else if inrange(`h_l',126,175) {
					local xlab = "150 200 250 300 350 400 450"
				}
				else {
					local xlab = ""
				}
				rdplot y_resid_sample tot_p if temp_in_reg==1, ///
					c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
					graph_options( ///
					title("`title'", color(black) size(large)) ///
					ytitle("`ytitle'", size(medlarge)) ///
					xtitle("2001 village population", size(medlarge)) ///
					ylabel(,nogrid angle(0) labsize(medlarge)) ///
					xlabel(`xlab', labsize(medlarge)) ///
					graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
					legend(off))
				graph export "$results/`folder'/`yvar'_fs_manual`h'_inreg_`bins'.pdf", replace
			}

		}	
	}
}


	// Append falsification results to VD elec first stage results file
keep fs_step-fdesc
dropmiss, obs force
tempfile bw_sens
save `bw_sens'
use "$results/RDROBUST_fs_vdelec.dta", clear
append using `bw_sens'
duplicates drop
compress
save "$results/RDROBUST_fs_vdelec.dta", replace

}

****************************************************************** 
****************************************************************** 

** 4. Placebo test for main VD_elec RDs
{
	//  Export data 
use "$panel/panel_dataset_full.dta", clear

	// Keep villages in RD sample
gen in_fs_sample = vplan4<11 & corr_state==1 & sample==1 & sample_h==1 & sing_h==1 & pop_non_zero==1 
keep if in_fs_sample==1
keep if pop_mismatch20==0 
	
	// Create state FEs (since RD robust doesn't let you pass them through)
drop if st_code==32 // Kerala, only 3 villages
tab st_code, gen(STFE)	
drop STFE1 // to avoid collinearity
tab st_code, gen(STFEH)	
drop STFEH1 STFEH2 STFEH12 // to avoid collinearity

	// Create lights-difference variable, to identify crazy outliers
gen lights_diff = abs(lights_max2011_hat - lights_max2001_hat)
keep if lights_diff<20

	// Rename variable to make them easier to loop over
rename vdp_pwr_d_com_11 vd_pwr_d_com_11	
rename vdp_pwr_h_*_avg_11 vd_pwr_h_*_11

	// Keep essential variables only
keep vd_pwr_d_com_11 vd_pwr_d_all_11 vd_pwr_h_com_11 tot_p lights_max1998 lights_max1999 lights_max2000 ///
	lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005 STFE* in_fs_sample

	// Save dataset for placebo tests
save "$sens/RDROBUST_placebo_test_vdelec.dta", replace
	
	// Run "PLACEBO_test_rdrobust_vdelec.do"

}

****************************************************************** 
****************************************************************** 






